library(KernSmooth)
library(xtable)
library(mgcv)
library(stargazer)
library(abind)

###############
## Load Data ##
###############

useCircuits <- 4:11
doAnalysis <- TRUE

round2 <- function(x) format(round(x,2),nsmall=2)
shadebands <- function(x,lo,hi,col) polygon(c(x,rev(x)),y=c(lo,rev(hi)),col=col,border=NA)

# load full judge dataset
judgedf <- read.delim("fjcdata.txt")
judgedf <- judgedf[,c("Judge.Identification.Number","Judge.Last.Name","Judge.First.Name","Judge.Middle.Name","Suffix","Birth.year","Gender","Race.or.Ethnicity","Court.Type","President.name","Party.Affiliation.of.President","Court.Type..2.","President.name..2.","Party.Affiliation.of.President..2.","Court.Type..3.","President.name..3.","Party.Affiliation.of.President..3.")]

# load chief judge dataset
chiefdf <- read.csv("chief_judge_party_by_circuit.csv")

# load surname data for race balance check
load("surnames2000.RData")

# create some descriptive stats across the full set of circuits
case_count_by_circuit <- rep(NA,11)
state_table_list <- list()
for (circuit in useCircuits){
    df <- read.csv(file=paste0("final_",circuit,"_df.csv"))
    case_count_by_circuit[circuit] <- nrow(df)
    state_table <- data.frame(Circuit=circuit,
                              State=names(table(df$state)),
                         Petitioners=as.numeric(table(df$state)),
                         Relief.Denials=table(df$state,df$panelvote == "Grant of Relief")[,1],
                         Executions=table(df$state,df$executed)[,2])
    rownames(state_table) <- NULL
    state_table_list[[which(circuit == useCircuits)]] <- state_table
}
state_table_all <- do.call(rbind,state_table_list)
print(xtable(state_table_all,digits=0),include.rownames=FALSE)

# create object to hold all models
data_objects <- list()
outcome_tables <- list()
fitted_models <- list()
balance_tests <- list()
balance_tests_race <- list()
treatment_year_tables <- list()
selection_effect_size_needed <- list()
panel_type_table <- matrix(NA,11,4,
                                     dimnames=list(1:11,c("DDD","DDR","DRR","RRR")))

for (circuit in useCircuits){
    # for testing # circuit <- 4
    
    df <- read.csv(file=paste0("final_",circuit,"_df.csv"))
    
    keepjudgeids <- unique(c(df$j1code,df$j2code,df$j3code))
    jdf <- judgedf[is.element(judgedf$Judge.Identification.Number,keepjudgeids),]
    
    # merge surname - race data
    
    surname_locs <- match(as.character(df$defendant.last),as.character(tolower(surnames2000$surname)))
    race_probs <- surnames2000[surname_locs,2:6]
    colnames(race_probs) <- c("white","black","hispanic","asian","other")
    
    # convert data to format for fitting models
    
    Y <- as.numeric(df$panelvote != "Grant of Relief")
    Y_lower <- as.numeric(df$lower.vote != "Grant of Relief")
    Z <- as.numeric(df$executed)
    M <- length(Y)
    Ij <- as.matrix(df[,c("j1code","j2code","j3code")])
    Statej <- df$state
    JudgeIDs <- sort(unique(c(df$j1code,df$j2code,df$j3code)))
    JudgeNames <- jdf$Judge.Last.Name[match(JudgeIDs,jdf$Judge.Identification.Number)]
    
    # Use USCA appointing pres party
    JudgeParty1 <- jdf$Party.Affiliation.of.President[match(JudgeIDs,jdf$Judge.Identification.Number)]
    JudgeParty2 <- jdf$Party.Affiliation.of.President..2.[match(JudgeIDs,jdf$Judge.Identification.Number)]
    JudgeParty <- replace(JudgeParty1,JudgeParty2 != "Reassignment" & jdf$Court.Type..2.[match(JudgeIDs,jdf$Judge.Identification.Number)] == "USCA",JudgeParty2[JudgeParty2 != "Reassignment" & jdf$Court.Type..2.[match(JudgeIDs,jdf$Judge.Identification.Number)] == "USCA"])
    
    JudgeCol <- c("blue","red")[1 + (JudgeParty == "Republican")]
    N <- length(JudgeIDs)
    Ij <- data.frame(factor(Ij[,1],levels= JudgeIDs),
                     factor(Ij[,2],levels= JudgeIDs),
                     factor(Ij[,3],levels= JudgeIDs))
    names(Ij) <- c("J1","J2","J3")
    JudgeCount <- tabulate(c(as.numeric(Ij[,1]),as.numeric(Ij[,2]),as.numeric(Ij[,3])))
    Treatment <- (JudgeParty[as.numeric(Ij[,1])] == "Democratic") + (JudgeParty[as.numeric(Ij[,2])] == "Democratic") + (JudgeParty[as.numeric(Ij[,3])] == "Democratic") < 2    
    CaseYear <- as.numeric(df$year)
    ExecutionYear <- as.numeric(df$executed.year)
    
    # construct chief judge party variable
    
    chiefdf_tmp <- chiefdf[,c(1,grep(circuit,colnames(chiefdf)))]
    colnames(chiefdf_tmp) <- c("year","chief_party")
    ChiefJudgeParty <- chiefdf_tmp$chief_party[match(CaseYear,chiefdf_tmp$year)]
    
    # descriptive stats
    
    tmptable <- table(Y,Z)
    dimnames(tmptable) <- NULL
    rownames(tmptable) <- c("Granted Relief","Denied Relief")
    colnames(tmptable) <- c("Not Executed","Executed")
    outcome_tables[[circuit]] <- tmptable
    print(tmptable)
    print(prop.table(tmptable,1))
    
    # diagnostic plot for outcomes over time
    
    pdf(paste0("OutcomesCaseYear",circuit,".pdf"),width=8,height=5)
    plot(CaseYear,Y,type="n",ylim=c(0,1),xlim=c(1977,2016),ylab="Proportion",xlab="Year",main=paste("Circuit",circuit))
    lines(ksmooth(CaseYear,Treatment,kernel="box",bandwidth=5),lwd=2,col="black")
    lines(ksmooth(CaseYear,Y,kernel="box",bandwidth=5),lwd=2,lty=2,col="grey")
    lines(ksmooth(CaseYear,Z,kernel="box",bandwidth=5),lwd=2,lty=3,col="grey")
    legend("topleft",legend=c("Rep Majority","Relief Denied","Executed"),lwd=c(2,2,2),lty=c(1,2,3),col=c("black","grey","grey"))
    dev.off()
    
    # create data object for model fitting
    
    df2 <- data.frame(
     Circuit = circuit,
     Denied = Y,
     Executed = Z,
     Republican_Majority = Treatment,
     Republican_Count = c("RRR","DRR","DDR","DDD")[(JudgeParty[as.numeric(Ij[,1])] == "Democratic") + (JudgeParty[as.numeric(Ij[,2])] == "Democratic") + (JudgeParty[as.numeric(Ij[,3])] == "Democratic") + 1],
     CaseYear = CaseYear,
     ChiefJudgeParty = ChiefJudgeParty,
     State = Statej,
     racep = race_probs,
     row.names = NULL
    )
    
    data_objects[[circuit]] <- df2
    
    panel_type_table[circuit,] <- table(df2$Republican_Count)
    treatment_year_tables[[circuit]] <- table(df2$Republican_Majority,df2$CaseYear)
    selection_effect_size_needed[[circuit]] <- diff(prop.table(table(df2$Denied,df2$Republican_Majority),1)[,2])
    
    # test for whether Republican_Majority predicts coming from a particular state
    
    for (st in unique(df2$State)){
      balance_tests[[st]] <- list(
      no_fe_p = summary(lm(I(State==st)~Republican_Majority,data=df2))$coefficients[2,4],
      fe_p = summary(lm(I(State==st)~Republican_Majority+as.factor(CaseYear),data=df2))$coefficients[2,4]
      )
    }
    
    # test for whether Republican_Majority is predicted by race|surname
    
    balance_tests_race[[circuit]] <- list(
    summary(lm(Republican_Majority~racep.black+as.factor(CaseYear),data=df2))$coefficients[2,4],
    summary(lm(Republican_Majority~racep.hispanic+as.factor(CaseYear),data=df2))$coefficients[2,4])
    
    # fit models
    
    model_denied_0 <- lm(Denied~Republican_Majority,data=df2)
    model_denied_fe <- lm(Denied~Republican_Majority+as.factor(CaseYear),data=df2)
    model_denied_2wayfe <- lm(Denied~Republican_Majority+as.factor(CaseYear)+as.factor(State),data=df2)
    
    model_executed_0 <- lm(Executed~Republican_Majority,data=df2)
    model_executed_fe <- lm(Executed~Republican_Majority+as.factor(CaseYear),data=df2)
    model_executed_2wayfe <- lm(Executed~Republican_Majority+as.factor(CaseYear)+as.factor(State),data=df2)
    
    model_denied_0_alt <- lm(Denied~relevel(Republican_Count,ref="DDR"),data=df2)
    model_denied_fe_alt <- lm(Denied~relevel(Republican_Count,ref="DDR")+as.factor(CaseYear),data=df2)
    model_denied_2wayfe_alt <- lm(Denied~relevel(Republican_Count,ref="DDR")+as.factor(CaseYear)+as.factor(State),data=df2)
    
    model_executed_0_alt <- lm(Executed~relevel(Republican_Count,ref="DDR"),data=df2)
    model_executed_fe_alt <- lm(Executed~relevel(Republican_Count,ref="DDR")+as.factor(CaseYear),data=df2)
    model_executed_2wayfe_alt <- lm(Executed~relevel(Republican_Count,ref="DDR")+as.factor(CaseYear)+as.factor(State),data=df2)
    
    if (circuit != 7){
    model_denied_0_split <- lm(Denied~Republican_Majority*ChiefJudgeParty,data=df2)
    model_denied_fe_split <- lm(Denied~Republican_Majority*ChiefJudgeParty+as.factor(CaseYear),data=df2)
    model_denied_2wayfe_split <- lm(Denied~Republican_Majority*ChiefJudgeParty+as.factor(CaseYear)+as.factor(State),data=df2)
    
    model_executed_0_split <- lm(Executed~Republican_Majority*ChiefJudgeParty,data=df2)
    model_executed_fe_split <- lm(Executed~Republican_Majority*ChiefJudgeParty+as.factor(CaseYear),data=df2)
    model_executed_2wayfe_split <- lm(Executed~Republican_Majority*ChiefJudgeParty+as.factor(CaseYear)+as.factor(State),data=df2)
    }
    
    # stargazer(model_0,model_fe,model_smooth_gaussian,model_smooth_logistic)
    
    if (circuit != 7){ fitted_models[[circuit]] <- list(model_denied_0,model_denied_fe,model_denied_2wayfe,
                                     model_executed_0,model_executed_fe,model_executed_2wayfe,
                                     model_denied_0_alt,model_denied_fe_alt,model_denied_2wayfe_alt,
                                     model_executed_0_alt,model_executed_fe_alt,model_executed_2wayfe_alt,model_denied_0_split,model_denied_fe_split,model_denied_2wayfe_split,
                                     model_executed_0_split,model_executed_fe_split,model_executed_2wayfe_split) } else {
        fitted_models[[circuit]] <- list(model_denied_0,model_denied_fe,model_denied_2wayfe,
                                     model_executed_0,model_executed_fe,model_executed_2wayfe,
                                     model_denied_0_alt,model_denied_fe_alt,model_denied_2wayfe_alt,
                                     model_executed_0_alt,model_executed_fe_alt,model_executed_2wayfe_alt)                               
                                     }
    

    
} # end loop over circuits

# Check balance tests

no_fe_p <- sapply(balance_tests,function(x) x$no_fe_p)
fe_p <- sapply(balance_tests,function(x) x$fe_p)
ks.test(no_fe_p,punif)
ks.test(fe_p,punif)

t.test(no_fe_p,mu=0.5)
t.test(fe_p,mu=0.5)

# Create coefficient tables for all circuits

coef_table <- se_table <- matrix(NA,max(useCircuits),6)
rownames(coef_table) <- rownames(se_table) <- 1:11
colnames(coef_table) <- colnames(se_table) <- c("Denial: Constant Time","Denial: Time FE","Denial: Time & State FE","Execution: Constant Time","Execution: Time FE","Execution: Time & State FE")
for (circuit in useCircuits){
  coef_table[circuit,1] <- coef(fitted_models[[circuit]][[1]])[2]
  coef_table[circuit,2] <- coef(fitted_models[[circuit]][[2]])[2]
  coef_table[circuit,3] <- coef(fitted_models[[circuit]][[3]])[2]
  coef_table[circuit,4] <- coef(fitted_models[[circuit]][[4]])[2]
  coef_table[circuit,5] <- coef(fitted_models[[circuit]][[5]])[2]
  coef_table[circuit,6] <- coef(fitted_models[[circuit]][[6]])[2]
}
se <- function(fit) sqrt(diag(vcov(fit)))
for (circuit in useCircuits){
  se_table[circuit,1] <- se(fitted_models[[circuit]][[1]])[2]
  se_table[circuit,2] <- se(fitted_models[[circuit]][[2]])[2]
  se_table[circuit,3] <- se(fitted_models[[circuit]][[3]])[2]
  se_table[circuit,4] <- se(fitted_models[[circuit]][[4]])[2]
  se_table[circuit,5] <- se(fitted_models[[circuit]][[5]])[2]
  se_table[circuit,6] <- se(fitted_models[[circuit]][[6]])[2]
}

# national estimates

national_total_estimates <- c(
  sum(coef_table[,1]*case_count_by_circuit,na.rm=TRUE),
  sum(coef_table[,2]*case_count_by_circuit,na.rm=TRUE),
  sum(coef_table[,3]*case_count_by_circuit,na.rm=TRUE),
  sum(coef_table[,4]*case_count_by_circuit,na.rm=TRUE),
  sum(coef_table[,5]*case_count_by_circuit,na.rm=TRUE),
  sum(coef_table[,6]*case_count_by_circuit,na.rm=TRUE)
)

national_proportion_estimates <- c(
  sum(coef_table[,1]*case_count_by_circuit,na.rm=TRUE)/sum(case_count_by_circuit,na.rm=TRUE),
  sum(coef_table[,2]*case_count_by_circuit,na.rm=TRUE)/sum(case_count_by_circuit,na.rm=TRUE),
  sum(coef_table[,3]*case_count_by_circuit,na.rm=TRUE)/sum(case_count_by_circuit,na.rm=TRUE),
  sum(coef_table[,4]*case_count_by_circuit,na.rm=TRUE)/sum(case_count_by_circuit,na.rm=TRUE),
  sum(coef_table[,5]*case_count_by_circuit,na.rm=TRUE)/sum(case_count_by_circuit,na.rm=TRUE),
  sum(coef_table[,6]*case_count_by_circuit,na.rm=TRUE)/sum(case_count_by_circuit,na.rm=TRUE)
)

se_sum <- function(x,na.rm=TRUE) sqrt(sum(x^2,na.rm=na.rm))
national_total_ses <- c(
  se_sum(se_table[,1]*case_count_by_circuit,na.rm=TRUE),
  se_sum(se_table[,2]*case_count_by_circuit,na.rm=TRUE),
  se_sum(se_table[,3]*case_count_by_circuit,na.rm=TRUE),
  se_sum(se_table[,4]*case_count_by_circuit,na.rm=TRUE),
  se_sum(se_table[,5]*case_count_by_circuit,na.rm=TRUE),
  se_sum(se_table[,6]*case_count_by_circuit,na.rm=TRUE)
)

national_proportion_ses <- c(
  se_sum(se_table[,1]*case_count_by_circuit,na.rm=TRUE)/sum(case_count_by_circuit,na.rm=TRUE),
  se_sum(se_table[,2]*case_count_by_circuit,na.rm=TRUE)/sum(case_count_by_circuit,na.rm=TRUE),
  se_sum(se_table[,3]*case_count_by_circuit,na.rm=TRUE)/sum(case_count_by_circuit,na.rm=TRUE),
  se_sum(se_table[,4]*case_count_by_circuit,na.rm=TRUE)/sum(case_count_by_circuit,na.rm=TRUE),
  se_sum(se_table[,5]*case_count_by_circuit,na.rm=TRUE)/sum(case_count_by_circuit,na.rm=TRUE),
  se_sum(se_table[,6]*case_count_by_circuit,na.rm=TRUE)/sum(case_count_by_circuit,na.rm=TRUE)
)

national_total_intervals <- rbind(national_total_estimates - 1.96*national_total_ses,national_total_estimates + 1.96*national_total_ses)
rownames(national_total_intervals) <- c("LB","UB")
colnames(national_total_intervals) <- colnames(coef_table)

# clean up for output

round2 <- function(x) format(round(x,2),nsmall=2)
addparens <- function(x) paste0("(",x,")")

tidy_coefs <- round2(coef_table[4:11,])
tidy_ses <- matrix(addparens(round2(se_table[4:11,])),8,6)
tidy_aggregate <- rbind(round2(national_proportion_estimates),addparens(round2(national_proportion_ses)),round(national_total_estimates),addparens(round(national_total_ses)),paste0(round(national_total_intervals[1,]),"-",round(national_total_intervals[2,])))
rownames(tidy_aggregate) <- c("Average Effect","(Standard Error)","Total Effect","(Standard Error)","95% Interval")

colnames(tidy_coefs) <- colnames(tidy_ses) <- colnames(tidy_aggregate) <- colnames(coef_table)
  
write.table(tidy_coefs,file=paste0("Coefs.tex"),sep="&",quote=FALSE)
write.table(tidy_ses,file=paste0("SEs.tex"),sep="&",quote=FALSE)
write.table(tidy_aggregate,file=paste0("Aggregates.tex"),sep="&",quote=FALSE)

# create coefficient plot

coef_table
se_table

pdf("CircuitEstimates.pdf",width=10,height=5)
par(mfrow=c(1,2))

plot(0,0,type="n",axes=FALSE,ylim=c(3.5,11.5),xlim=c(-0.4,0.4),xlab="",ylab="",main="Denial")
axis(1)
axis(2,at=4:11,las=2)
rect(national_proportion_estimates[2] + -1.96*national_proportion_ses[2],
     3.5,
     national_proportion_estimates[2] + 1.96*national_proportion_ses[2],
     11.5,
     col=rgb(0,0,0,0.2),border=NA)
lines(rep(national_proportion_estimates[2],2),c(3.5,11.5),col=rgb(0,0,0,0.5),lwd=2)
for (i in 4:11) points(coef_table[i,2],i,pch=16)
for (i in 4:11) lines(coef_table[i,2]+c(-1.96,1.96)*se_table[i,2],c(i,i),pch=16)

plot(0,0,type="n",axes=FALSE,ylim=c(3.5,11.5),xlim=c(-0.4,0.4),xlab="",ylab="",main="Execution")
axis(1)
axis(2,at=4:11,las=2)
rect(national_proportion_estimates[5] + -1.96*national_proportion_ses[5],
     3.5,
     national_proportion_estimates[5] + 1.96*national_proportion_ses[5],
     11.5,
     col=rgb(0,0,0,0.2),border=NA)
lines(rep(national_proportion_estimates[5],2),c(3.5,11.5),col=rgb(0,0,0,0.5),lwd=2)
for (i in 4:11) points(coef_table[i,5],i,pch=16)
for (i in 4:11) lines(coef_table[i,5]+c(-1.96,1.96)*se_table[i,5],c(i,i),pch=16)
dev.off()

# what proportion of cases are in a circuit-year where there are no D majority or no R majority cases?

tidy_fe_informative <- t(t(round2(unlist(lapply(treatment_year_tables,function(x) sum(x[,(x[1,] != 0) & (x[2,] != 0)])/sum(x))))[4:11]))
rownames(tidy_fe_informative) <- rownames(tidy_coefs)

write.table(tidy_fe_informative,file=paste0("FE_informative.tex"),sep="&",quote=FALSE)

tidy_fe_informative <- sum(unlist(lapply(treatment_year_tables,function(x) sum(x[,(x[1,] != 0) & (x[2,] != 0)])))) / sum(unlist(lapply(treatment_year_tables,function(x) sum(x[,]))))

# Create coefficient tables for all circuits for alternative analysis with four treatment levels

coef_table <- se_table <- array(0,c(max(useCircuits),6,4))
dimnames(coef_table)[[1]] <- dimnames(se_table)[[1]] <- 1:11
dimnames(coef_table)[[2]]  <- dimnames(coef_table)[[2]] <- c("Denial: Constant Time","Denial: Time FE","Denial: Time & State FE","Execution: Constant Time","Execution: Time FE","Execution: Time & State FE")
dimnames(coef_table)[[3]] <- dimnames(se_table)[[3]] <- c("DDD","DDR","DRR","RRR")
for (circuit in useCircuits){
  coef_table[circuit,1,c(1,3,4)] <- coef(fitted_models[[circuit]][[7]])[2:4]
  coef_table[circuit,2,c(1,3,4)] <- coef(fitted_models[[circuit]][[8]])[2:4]
  coef_table[circuit,3,c(1,3,4)] <- coef(fitted_models[[circuit]][[9]])[2:4]
  coef_table[circuit,4,c(1,3,4)] <- coef(fitted_models[[circuit]][[10]])[2:4]
  coef_table[circuit,5,c(1,3,4)] <- coef(fitted_models[[circuit]][[11]])[2:4]
  coef_table[circuit,6,c(1,3,4)] <- coef(fitted_models[[circuit]][[12]])[2:4]
}
coef_table[1:3,,] <- NA
for (circuit in useCircuits){
  se_table[circuit,1,c(1,3,4)] <- se(fitted_models[[circuit]][[7]])[2:4]
  se_table[circuit,2,c(1,3,4)] <- se(fitted_models[[circuit]][[8]])[2:4]
  se_table[circuit,3,c(1,3,4)] <- se(fitted_models[[circuit]][[9]])[2:4]
  se_table[circuit,4,c(1,3,4)] <- se(fitted_models[[circuit]][[10]])[2:4]
  se_table[circuit,5,c(1,3,4)] <- se(fitted_models[[circuit]][[11]])[2:4]
  se_table[circuit,6,c(1,3,4)] <- se(fitted_models[[circuit]][[12]])[2:4]
}
se_table[1:3,,] <- NA

# translate to national-level estimates
national_total_estimates <- array(NA,dim(coef_table)[2:3])
for (k in 1:4) national_total_estimates[,k] <- c(
  sum(coef_table[,1,k]*case_count_by_circuit,na.rm=TRUE),
  sum(coef_table[,2,k]*case_count_by_circuit,na.rm=TRUE),
  sum(coef_table[,3,k]*case_count_by_circuit,na.rm=TRUE),
  sum(coef_table[,4,k]*case_count_by_circuit,na.rm=TRUE),
  sum(coef_table[,5,k]*case_count_by_circuit,na.rm=TRUE),
  sum(coef_table[,6,k]*case_count_by_circuit,na.rm=TRUE)
)
national_proportion_estimates <- array(NA,dim(coef_table)[2:3])
for (k in 1:4) national_proportion_estimates[,k] <- c(
  sum(coef_table[,1,k]*case_count_by_circuit,na.rm=TRUE)/sum(case_count_by_circuit,na.rm=TRUE),
  sum(coef_table[,2,k]*case_count_by_circuit,na.rm=TRUE)/sum(case_count_by_circuit,na.rm=TRUE),
  sum(coef_table[,3,k]*case_count_by_circuit,na.rm=TRUE)/sum(case_count_by_circuit,na.rm=TRUE),
  sum(coef_table[,4,k]*case_count_by_circuit,na.rm=TRUE)/sum(case_count_by_circuit,na.rm=TRUE),
  sum(coef_table[,5,k]*case_count_by_circuit,na.rm=TRUE)/sum(case_count_by_circuit,na.rm=TRUE),
  sum(coef_table[,6,k]*case_count_by_circuit,na.rm=TRUE)/sum(case_count_by_circuit,na.rm=TRUE)
)
national_total_ses <- array(NA,dim(coef_table)[2:3])
for (k in 1:4) national_total_ses[,k] <- c(
  se_sum(se_table[,1,k]*case_count_by_circuit,na.rm=TRUE),
  se_sum(se_table[,2,k]*case_count_by_circuit,na.rm=TRUE),
  se_sum(se_table[,3,k]*case_count_by_circuit,na.rm=TRUE),
  se_sum(se_table[,4,k]*case_count_by_circuit,na.rm=TRUE),
  se_sum(se_table[,5,k]*case_count_by_circuit,na.rm=TRUE),
  se_sum(se_table[,6,k]*case_count_by_circuit,na.rm=TRUE)
)
national_proportion_ses <- array(NA,dim(coef_table)[2:3])
for (k in 1:4) national_proportion_ses[,k] <- c(
  se_sum(se_table[,1,k]*case_count_by_circuit,na.rm=TRUE)/sum(case_count_by_circuit,na.rm=TRUE),
  se_sum(se_table[,2,k]*case_count_by_circuit,na.rm=TRUE)/sum(case_count_by_circuit,na.rm=TRUE),
  se_sum(se_table[,3,k]*case_count_by_circuit,na.rm=TRUE)/sum(case_count_by_circuit,na.rm=TRUE),
  se_sum(se_table[,4,k]*case_count_by_circuit,na.rm=TRUE)/sum(case_count_by_circuit,na.rm=TRUE),
  se_sum(se_table[,5,k]*case_count_by_circuit,na.rm=TRUE)/sum(case_count_by_circuit,na.rm=TRUE),
  se_sum(se_table[,6,k]*case_count_by_circuit,na.rm=TRUE)/sum(case_count_by_circuit,na.rm=TRUE)
)


national_total_intervals <- abind(national_total_estimates - 1.96*national_total_ses,national_total_estimates + 1.96*national_total_ses,along=3)
dimnames(national_total_intervals)[[3]] <- c("LB","UB")
dimnames(national_total_intervals)[[1]] <- dimnames(coef_table)[[2]]
dimnames(national_total_intervals)[[2]] <- dimnames(coef_table)[[3]]


# clean up for output

tidy_aggregate <- paste(round2(national_proportion_estimates),apply(round2(national_proportion_ses),c(1,2),addparens))
tidy_aggregate <- matrix(tidy_aggregate,nrow=nrow(national_proportion_estimates),ncol=ncol(national_proportion_estimates))
rownames(tidy_aggregate) <- dimnames(coef_table)[[2]]
colnames(tidy_aggregate) <- dimnames(coef_table)[[3]]

write.table(panel_type_table,file=paste0("PanelTypeTable.tex"),sep="&",quote=FALSE)

write.table(tidy_aggregate,file=paste0("Aggregates4Level.tex"),sep="&",quote=FALSE)



# how big a selection effect would be needed?

selection_effect_required <- sum(unlist(selection_effect_size_needed)*case_count_by_circuit[4:11])/sum(case_count_by_circuit[4:11])




# Create coefficient tables for all circuits for alternative analysis with chief judge party interaction

chief_judge_interaction <- data.frame(circuit = 1:11,denial_est=NA,denial_se=NA,execution_est=NA,execution_se=NA)
for (circuit in setdiff(useCircuits,7)){
  
denial_coef_tmp <- coef(fitted_models[[circuit]][[13]])
denial_se_tmp <- se(fitted_models[[circuit]][[13]])
chief_judge_interaction[circuit,2:3] <- c(
                     denial_coef_tmp[names(denial_coef_tmp) == "Republican_MajorityTRUE:ChiefJudgePartyR"],
denial_se_tmp[names(denial_se_tmp) == "Republican_MajorityTRUE:ChiefJudgePartyR"])

execution_coef_tmp <- coef(fitted_models[[circuit]][[16]])
execution_se_tmp <- se(fitted_models[[circuit]][[16]])
chief_judge_interaction[circuit,4:5] <- c(
                     execution_coef_tmp[names(execution_coef_tmp) == "Republican_MajorityTRUE:ChiefJudgePartyR"],
execution_se_tmp[names(execution_se_tmp) == "Republican_MajorityTRUE:ChiefJudgePartyR"])
}

sum(case_count_by_circuit[setdiff(useCircuits,7)] * chief_judge_interaction$denial_est[setdiff(useCircuits,7)])/sum(case_count_by_circuit[setdiff(useCircuits,7)])

se_sum(chief_judge_interaction$denial_se[setdiff(useCircuits,7)]*case_count_by_circuit[setdiff(useCircuits,7)],na.rm=TRUE)/sum(case_count_by_circuit[setdiff(useCircuits,7)],na.rm=TRUE)

sum(case_count_by_circuit[setdiff(useCircuits,7)] * chief_judge_interaction$execution_est[setdiff(useCircuits,7)])/sum(case_count_by_circuit[setdiff(useCircuits,7)])

se_sum(chief_judge_interaction$execution_se[setdiff(useCircuits,7)]*case_count_by_circuit[setdiff(useCircuits,7)],na.rm=TRUE)/sum(case_count_by_circuit[setdiff(useCircuits,7)],na.rm=TRUE)




# calculate aggregate outcomes

aggregate_outcome_table <- outcome_tables[[11]] # create object of correct shape to fill in below
aggregate_outcome_table[1,1] <- sum(unlist(lapply(outcome_tables,function(x) x[1,1])))
aggregate_outcome_table[1,2] <- sum(unlist(lapply(outcome_tables,function(x) x[1,2])))
aggregate_outcome_table[2,1] <- sum(unlist(lapply(outcome_tables,function(x) x[2,1])))
aggregate_outcome_table[2,2] <- sum(unlist(lapply(outcome_tables,function(x) x[2,2])))

prop.table(colSums(aggregate_outcome_table))
prop.table(rowSums(aggregate_outcome_table))


# balance checks for race|surname

# by circuit
mean(unlist(lapply(balance_tests_race,function(x) x[[1]]))) # avg p-value for predicting treatment with black|surname prob across circuit
mean(unlist(lapply(balance_tests_race,function(x) x[[2]]))) # avg p-value for predicting treatment with hispanic|surname prob across circuit

# across all circuits
merged_df <- do.call(rbind,data_objects)

# balance check
summary(lm(Republican_Majority~racep.hispanic+racep.black+as.factor(CaseYear)*as.factor(Circuit),data=merged_df))$coefficients[1:3,]

# check that overall results do not change if we control for defendant race|surname

summary(lm(Denied~Republican_Majority+as.factor(CaseYear)*as.factor(Circuit),data=merged_df))$coefficients[1:2,]
summary(lm(Executed~Republican_Majority+as.factor(CaseYear)*as.factor(Circuit),data=merged_df))$coefficients[1:2,]

summary(lm(Denied~Republican_Majority+racep.hispanic+racep.black+as.factor(CaseYear)*as.factor(Circuit),data=merged_df))$coefficients[1:4,]
summary(lm(Executed~Republican_Majority+racep.hispanic+racep.black+as.factor(CaseYear)*as.factor(Circuit),data=merged_df))$coefficients[1:4,]




